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ABSTRACT 

We present an efficient method to compute CMB anisotropics in non-flat 
universes. First we derive the Boltzmann equation for cosmic microwave 
background temperature and polarization fluctuations produced by scalar 
perturbations in a general Robertson- Walker universe. We then apply the 
integral method to solve this equation, writing temperature and polarization 
anisotropics as a time integral over a geometrical term and a source term. The 
geometrical terms can be written using ultra-spherical Bessel functions, which 
depend on curvature. These cannot be precomputed in advance as in flat space. 
Instead we solve directly their differential equation for selected values of the 
multipoles. The resulting computational time is comparable to the flat space 
case and improves over previous methods by 2-3 orders of magnitude. This 
allows one to compute highly accurate CMB temperature and polarization 
spectra, matter transfer functions and their CMB normalizations for any 
cosmological model, thereby avoiding the need to use various inexact fitting 
formulae that exist in the literature. 



Subject headings: cosmology: cosmic microwave background, cosmology: 
large-scale structure of the universe, gravitation, cosmology: dark matter 
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1. Introduction 



In the past few years the field of cosmic microwave background (CMB) anisotropics 
transformed from a theoretical exercise to an active experimental area of research. Since 
the first discovery by COBE in 1992 (Smoot et al. 1992) there have been more than a 
dozen independent detections of fiuctuations over a much larger angular range (see e.g. 
Lineweaver et al. 1996 and Rocha and Hancock 1996 for compilations). The future looks 
even more promising, as there are two funded satellite missions in planning, one in Europe 
(Planck) and one in USA (MAP), that have the promise to measure the anisotropics to 
an exquisite accuracy over three decades in angular scale. This will open the possibility 
of determining many cosmological parameters with much greater accuracy than any other 
probe in cosmology (|Jungman et al. 1996 ; Bond, Efstatiou fc Tegmark 1997 ; Zaldarriaga, 
tSpergel fc Seljak 199q ). 

In light of this experimental promise theoretical predictions have been advanced to 
a higher stage as well. It is not sufficient anymore to make qualitative statements and 
approximate predictions to analyze the forthcoming data and to guide experimental design. 
Rather what is needed are very accurate predictions of CMB spectra (typically to better 
than 1%), which allow one to study the CMB sensitivity to various cosmological parameters. 
Because this sensitivity depends on the choice of mission specifics such as number of 
detectors, their noise characteristics or angular resolution, it is important that such accurate 
predictions are available already at the mission design stage when its characteristics are 
still open to modification. In Seljak & Zaldarriaga 1996 (hereafter Paper I) we presented 
a method that is both fast, accurate and applicable to any fiat cosmological model. The 
method is based on the source time integration over the photon past light cone and has the 
advantage of reducing the computational time by about two orders of magnitude compared 
to the more traditional methods, while still being exact in the sense that it can achieve 
arbitrary precision within the limits of linear perturbation theory. Such a method therefore 
allows one to explore a large range of parameter space with high accuracy. 

Open cosmological models and their predictions for CMB have received much attention 
lately (e.g. Gouda & Sugiyama 1992, Kamionkowski, Spergel & Sugiyama 1994, Hu, Bunn 
& Sugiyama 1995, White & Bunn 1995, Gorski et al. 1995). The reason for this is simple: 
observational evidence from the nearby universe (e.g. Peebles 1993) suggests that an open 
universe is favored over the critical or closed one. While the differences might be resolved 
with the addition of a cosmological constant A, so that the geometry of the universe would 
remain flat, the value of A is already significantly constrained by a number of independent 
tests based on COBE data, lensing statistics and SN type la results (Bunn & White 1997; 
Kochanek 1996; Perlmutter et al. 1996). It therefore seems natural to explore the possibility 
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that the universe is indeed open. From a theoretical point of view, open inflationary models 
proved to be constructible and capable of generating the initial perturbations (e.g. Lyth & 
Stewart 1990, Ratra & Peebles 1994, Liddle et al. 1995, Bucher et al. 1995). Models with 
positive curvature, although currently less popular, have also been studied in detail (White 
& Scott 1996). 

All the computations of CMB anisotropies in open and closed universes carried out so 
far have used the traditional Boltzmann hierarchy approach pioneered by Peebles & Yu 
(1970) and extended by Wilson & Silk (1981) and Bond & Efstathiou (1984). The speed 
limitations of the Boltzmann approach in fiat space were extensively discussed in Paper I. 
In the case of non-flat geometry they become even worse: solving the hierarchy in an open 
universe is much slower than in the fiat case and leads to extremely long integration times. 
This malces an accurate search over a large set of parameter space practically impossible. 
Instead one is forced to use approximations, but as shown in this paper these give at best a 
few percent accuracy and can be very misleading in the study of parameter determination. 

In this paper we develop the integral solution for photon transport in a general 
Robertson- Walker background, thereby generalizing the method developed in Paper I to a 
non-fiat geometry. As in the fiat case the solution is written in terms of a time integral over a 
source term and a geometrical term. The latter can be expressed in terms of ultra-spherical 
Bessel functions, defined as the radial part of the eigenfunction of the Laplacian on a curved 
manifold. We present a method for computing these functions efficiently in the context of 
CMB calculations and incorporate them into the integral solution. This allows us to achieve 
a fast and accurate method, improving over previous calculational methods by 2-3 orders of 
magnitude. In this paper we concentrate on the temperature and polarization anisotropies 
produced by scalar modes, the generalization to vector and tensor perturbations will be 
presented elsewhere ( [Hu, Seljak, White fc Zaldarriaga 1997| ). The outline of the paper is as 



follows: in §2 we present the Einstein and fluid equations in a general Robertson Walker 
background, in §3 we discuss the Boltzmann hierarchy for CMB photons and in §4 we derive 
its integral solution. In §5 we discuss in detail the method for computing the ultra-spherical 
Bessel functions. The comparison between the exact solution and the approximations often 
used in the literature is presented in §6. In §7 we discuss the numerical implementation of 
the method. This is followed by discussion and conclusions in §8. In an appendix we present 
an alternative derivation of the integral solution, highlighting its geometrical interpretation. 



2. Einstein and Fluid Equations 
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In this section we present the Einstein and fluid differential equations for the metric, 
cold dark matter (CDM) and baryons that must be solved to calculate the CMB anisotropy 
spectra. These equations are the basis of the traditional methods and are also used in the 
integral method, discussed in §4. The derivation of the Einstein and fluid equations in a 
non-flat background can be found in the literature (see for example Bertschinger 1996), so 
we just present the final results. 

The metric is written as 

ds^ = —dt^ + a^(7ij + hij)dx^dx^ 

= a'^l-dr'^ + {'^ij + hij)dx'dx^], (1) 

where a is the expansion factor, Xj the co moving coordinates and r = J dt/a the conformal 
time. We are using units in which c = 1. Space part of unperturbed metric is jij with 
constant curvature K = Hq{Qq — 1) and hij is the metric perturbation in synchronous 
gauge (Lifshitz 1946). Although all observable quantities are identical in different gauges 
the computational efficiency to obtain them within a given accuracy is not. This criterion 
lead us to work in synchronous gauge [|. In comparison to the longitudinal gauge (Pardeen 



1980| ) it is about 20% more efficient with isentropic initial conditions and even more so with 



isocurvature initial conditions, which are difficult to set up in the longitudinal gauge. 

When working with linear theory in a flat universe it is convenient to use Fourier 
modes as they evolve independently. Their generalization in a non-flat universe are the 
eigenfunctions of the Laplacian operator that we shall call G{k, x) (e.g. Abbott & Schaefer 
1986), 

\/^Gik,x) = -k^Gik,x). (2) 

We expand all the perturbations in terms of G and its spatial covariant derivatives. 
For example, the metric perturbations for a single mode are given by 

hij = ^7.,G +ih + 6ri){k-^G^,, + ^^,,G), (3) 

where h and {h + 6r]) are the trace and traceless part of the metric perturbation. Covariant 
derivatives with respect to the metric 7jj are denoted with a "|" The perturbed Einstein's 
equations result in the following equations for h and t] (Bertschinger 1996), 

{k^ - 3K)ri - l--h = -SnQa^Sp 
Zi a 

{k'^ - 3K)f] + = A7iga\p + p)kv, (4) 



In paper I we presented the equations in the longitudinal gauge. 
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Q here stands for the gravitational constant; Sp and v characterize the density and velocity 
perturbations [v = ik ■ v), 6p = J2j Pj^j-, (p + P)'v = HjiPj + Pj)'Vj-, where pj and pj are 
the mean density and pressure of the j-th species and the sum is carried out over all the 
different species in the universe. 

The equation for the cold dark matter density perturbation 5c is, 

4 = (5) 

where by definition in this gauge the cold dark matter particles have zero peculiar velocities. 
The Euler equation for the baryons has additional terms caused by Thomson scattering and 
pressure, so baryons have velocities relative to the dark matter, 

5b = -kvb - ^ , 

Vb = --Vb + clk5b + ^arieXeariv^ - Vb) . (6) 
a Spb 

Here Cg is the baryon sound speed, Vb is the baryon velocity, v^ is given by the temperature 
dipole v^ = 3Ati and p^, pb are the mean photon and baryon densities, respectively. The 
Thomson scattering cross section is ctt, rie is the electron density and Xe is the ionization 
fraction. 



3. Boltzmann equation 



In this section we discuss CMB anisotropics and present a derivation of the Boltzmann 
equation for the photons. Because we use a novel approach to treat this problem (see also 
Hu & White 1997) we start by considering flat geometries in §3.1 and then generalize 
the results to arbitrary Robert son- Walker backgrounds in §3.2. The derivation of the 
Boltzmann hierarchy for polarization is new to this work. The notation we use is largely 
based on Paper I. 



The CMB radiation field is described by a 2 x 2 intensity tensor lij ( Phandrasekar 



1960|) . The Stokes parameters Q and U are defined as Q = (In — /22)/4 and U = /12/2, 
while the temperature anisotropy is given by T = (In + /22)/4. The fourth Stokes 
parameter V that describes circular polarization is not necessary in standard cosmological 
models because it cannot be generated through the process of Thomson scattering. While 
the temperature is a scalar quantity Q and U are not. They depend on the direction of 
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observation n and on the two axis (ei, 62) perpendicular to n used to define them. If for a 
given h the axes (ei, 62) are rotated by an angle ■?/' such that e'^ = cos?/' ei + sin-?/; 62 and 
62 = — sin ip ei + cos 62 the Stokes parameters change as 

Q' = cos 2iIj Q + sin 2^/- f/ 

U' = - sin 27/- Q + cos 2?/' U (7) 



To analize the CMB temperature on the sky it is natural to expand it in spherical 
harmonics. These are not appropriate for polarization, because the two combinations 
Q ± iU are quantities of spin ±2 ( [Goldberg et al. 1967| ). They should be expanded in 
spin- weighted harmonics ±2^"^ ( jZaldarriaga fc Seljak 1997| ), 



Tin) = J2 o-TM^imin) 

Im 

{Q + iU){h) = '^a2,im 2Yim{n) 



Im 

{Q-iU){h) = '^a^2,im -^yimin). (8) 

Im 

There is an equivalent expansion using tensors on the sphere (Kamionkowski, Kosowsky & 
Stebbins 1997). The coefficients a±2,im are observable on the sky and their power spectra 
can be predicted for different cosmological models. Instead of 0-1-2,;™ it is convenient to use 
their linear combinations aE,im = -{ci2,im + a_2,«m)/2 and as^im = -{a2,im - a_2,«m)/2i, 
which have opposite parities. Four power spectra are needed to characterize fluctuations in 
a gaussian theory, the autocorrelation of T, E and B and the cross correlation of E and T. 
Because of parity considerations the cross-correlations between B and the other quantities 
vanish and one is left with 



Cxi 

Cci 



1 



21 + 
1 



21 + 



(9) 



where X stands for T, E or B and (■ ■ ■) means ensemble average. Only vector and tensor 
modes contribute to B ( [Zaldarriaga fc Seljak 1997| , [Kamionkowski et al. 199"7|) , hence we 
may ignore it in the remainder of this paper. 



3.1. Flat geometry 
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We will start by studying perturbations in a flat universe to introduce our notation and 
clarify the treatment of polarization, which differs from the usual method where Legendre 
polynomials are used to expand both temperature and polarization (e.g. Bond & Efstathiou 
1984). As we are dealing with a linear problem we may consider only one eigenmode of 
the Laplacian (a plane wave in the fiat case) at a time. We may choose without loss of 
generality that || z. To define the Stokes parameters we use the spherical coordinate unit 
vectors {eg,e^). In this particular coordinate system only Q is different from zero and we 
denote it by Ap, so that Ap = Q = Q ± ill {U = 0). The temperature anisotropy for the 
single eigenmode is denoted by A^. 

For a single plane wave rotational symmetry implies that both Ay and Ap depend 
only on the angle between n and z {k || z), so only harmonics with m = are needed in the 
expansion. To calculate the evolution of these two quantities, we expand them as 



AT{k,n) 



A, 



G{k,x)J2{ 
Gik,x)J2{ 

G(fc,f)E( 
G(fc,f)E(- 



—I 



—I 



—I 



47r(2/ + 1)Ati 
(2/ + l)ATiP/(/i) 



47r(2/ + l)(/ + 2)!/(/-2)!2A„ ^Y^i^i) 



47r(2/ + !)(/ + 2)!/(/ - 2)! _2A„ -s^A/i) 
(2/ + 1) 2A„Pf(^) 



(10) 



where G{k, x) = exp{ik ■ x) and fi = k ■ n. We added a subindex ±2 to ±2^Pi to denote 
that they are the expansion coefficients in spin ±2 harmonics and we used the explicit 
expression for spin s harmonics with m = to write them in terms of associated Legendre 
polynomials ( [Goldberg et al. 1967D , 



\2l + l] 



Pi{cos9) 



±2^1 



(2/ + 1) (/-2)! 
\ 47r (/ + 2)! 



P/Ccosi 



As stated above for scalar modes in this reference frame U = 0, so Ap describes both 
spin ±2 quantities. For m = one has 2Y^ = -2^° and so 2Ap; = _2^p/- For density 



^The relation between these coefficient and those used in Zaldarriaga & Seljak (1997) is ±2Api = 
-^{l^2)l/il + 2)lAEi. 
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perturbations only m = harmonics are needed in the expansion because of the azimuthal 
symmetry, but the treatment can be generahzed to vector and tensor modes. In these cases 
modes with m = ±1 and m = ±2 are required, respectively. For vector and tensor modes 
U is no longer zero, so separate expansions for both Q + ill and Q — ill are needed. 

The Boltzmann equation for the CMB photons reads (e.g. |Ma fc Bertschinger 1995| ), 



At + ikfiAx 
Ap + ikfiAp 



o o 

A 



PlThomson- 



(12) 



The first term in the temperature equation represents the effect of gravitational redshift, 
while Arp\Thomson ^ud Ap^Thomson ^rc the changes in the photon distribution function 
produced by Thomson scattering, where the derivatives are taken with respect to conformal 
time. After inserting equation (|10D into equation ([l^) one obtains a system of two coupled 
hierarchies, 



Ato 
Ati 

At2 

Ati 

2Api 



-kA 



Tl 



h . 

— + ATOlThomson 
D 



[A 



TO 



2A 



T2\ 



2A. 



is) 

Tl 



3A 



T3 



ATl\Thomson 

2 2 • 

+ —k a + AT2\Thomson 
15 



k 



21 + 1 

k 



IAt(i-i) — {I + 1)At{i+i) 



+ A 



Tl\Thomson 



l>2 



[(/ - 2) 2Api_i - (/ + 3) 2A„+i] + 2Api\Thc 



21 + 1 

where a = {h + 6i))/2k'^, and we used the recurrence relations for the Legendre functions 

1 



(13) 



21 + 



-[/Pz_i + (/ + 1)Phi 



21+ 1 



(/ + 2)Pf_i + (/ - l)p2^i 



(14) 



The Thomson scattering cross section is. 



da 

dn 



e ■ e 



in 



(15) 



where e and e' are the unit vectors that describe the polarization of the electric field of the 
scattered and incoming radiation, respectively. The scattering terms in equations ( [T3| ) are 
most easily computed in the coordinate system where the incident photons travel along 
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the z axis and the electrons are at rest. If n' is the direction of the incident photon and h 
that of the scattered one then n' = i = = 0, = 0) and (6', 0) describe n. For a given 
scattering event, the Thomson scattering matrix is the simplest when expressed in terms 
of the intensities of radiation parallel (Ty) and perpendicular (Ti) to the plane containing 
both h and f\! . Equation ([T5|) leads to the following relation between incoming and scattered 
radiation, 



f|| = — arcos^^fi; 
Stt " 

Svr 



(16) 



where gt is the Thomson scattering cross section. The total intensity is the sum of the two 
components, T = T|| + T^, while the difference gives polarization Q = Ty — T^. Because 
the components are measured using this coordinate system the Stokes parameters of the 
incoming radiation Q' and U' depend on the angle of the scattered photon, while Q 
and V are already measured relative to the correct frame. It is more useful to refer the 
Stokes parameters of the incoming radiation relative to a fixed frame. To achieve this 
we construct the scattering matrix in terms of T', Q' + iV = exp{2i(j)){Q' + ill') and 
Q' — iU' = exp{—2i(j)){Q' — iU'), where we have used the transformation law (equation |^) 
to relate the two sets of Stokes parameters. 

Equation ([T6| ) implies that the scattered radiation in direction h is 

-- ^ [7(1 + cos^ e)T' + ^(cos^ e - l)e'''^{Q' + tU') + 



6T{h', n) 



5{Q±zU){n',n) = ^ 

47r 



-(cos^ e - 1)T' + -{cose± lfe-^'^{Q' + iU') + 
4 8 



8 



The final expression for the scattered field is an integral over all directions n'. 



(17) 



Xifi) 



Thomson 



—aaxne Xe 



X{n) + / dn'6X{h\h] 



where X stands for T and {Q ± ill). The first term accounts for the photons that are 
scattered away from the line of sight and the expansion factor a is introduced because we 
are calculating the derivative with respect to conformal time. 
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Equation ([T7D for the scattering matrix is written in the frame where 



n 



0, 0' = 0). We can use equation (|TT]) to show that ±2^2^ 



n 



^(lTcos^)V2e 



2i4 



and ±2^2 ^( 



n] 



An 



1 ± COS ( 



2^2e _ These together with the exphcit expressions 
^ 5mn. and +2Y^(n') = J T-Sm:f2 unable us to rewrite (|l 



in a more useful form {6ij is the Kronecker delta), 



6T{n, h) 



^2Yr{n'){Q' -lU') 



5{Q±iU){n\h) 



^- ^.Y^in) oYr{n)r + ^ ^.Y^ih) 2Y,-{n'){Q' + zf/') + 



-2-^2 



;n')(Q'-z^': 



(19) 



The advantage of this form for the scattering matrix comes from the fact that we want the 



scattering matrix in the frame where k \\ z and not h' = z. The sum J2m sYi"^{h) 



s'J-l 



n') 



acquires a phase change under rotation of the coordinate system that exactly cancels the 
phase change in the transformation of {Q ± iU) in equation (|19D, we may therefore use this 
equation in the coordinate system where A; || z to compute the Thomson scattering terms in 
equation ([T3|). Equation (|1^) is also particularly useful to perform the integral in equation 



Substituting the expansion for the Stokes parameters from equation (|TD|) into equation 
(P!P|) and using equation (|T^) we find 



A. 



Tl \Thomson 



Ati + J dnoYf'{n)ST{n) 



n 



±2^Pl\Thomson = —aarXeUe 



±2A« + j dn^Yrih) 6{Q±tU){n) 



n 



k{ ±2Api - —S12) 



(20) 



with n = At2 — 6( 2'Ap; + _2Ap;). The differential optical depth for Thomson scattering 
is denoted k = aUeXeCT- Note that the polarization has sources only at / = 2. One can 
apply the same analysis to vector and tensor perturbations, the only difference is that in 
equation (|TUp the expansion is in terms of harmonics with m = ±1, ±2 for vectors and 
tensors, respectively. In fact, the form of the scattering terms in equation (pO]) is the same 
for the three types of perturbations, as also noted by Hu & White (1997). Equation (120) is 
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valid in the rest frame of the electrons, so in the reference frame where the baryon velocity 
is ff, the distribution of scattered radiation has an additional dipole. The final expression 
for the Boltzmann hierarchy is 



A TO 

Ati 
At; 

2^ PI 



-kA 



Tl 



3 [Ato 



2At21 +k[j-ATi 



2,A^-^ — 3A'P3 



Vb 

3 

2 ,2 

— k a + K 
15 



n 

10 



k 



21 + 1 
k 



IAt{i-i) — {I + 1)At{i+i) 



kAti J > 2 



21 



[{I - 2) 2A„_i - (/ + 3) 2A„+i] 



k ^Api - —kliSn. (21) 



3.2. Non-fiat geometry 



We now proceed to generalize the results of the previous section to a general 
Robert son- Walker background. First we generalize equation (^). Following Wilson & Silk 
(1981) we expand the photon temperature At in terms of Legendre tensors, 

At = E(2/ + 1) ATi{-Pr\\{h)-'G\,,...,, Pr\ (22) 
I I 

with 13"^ = k'^ + K , h'f = 1 — / fS"^ and \ii ■ ■ - ii covariant derivatives. These Legendre 
tensors are symmetric combinations constructed out of n and the metric g"^^ using the 
recursion 

{21 + l)n^' Pi'-''^ = {l + l)PiXi'^' + Ig^''' Pt'^''\ (23) 

where the parentheses imply symmetrization with respect to the indices (Wilson & Silk 
1981; [White fc Scott 1995|) . The first moment at / = is given by Pq = 1- In the flat 



universe this recursion reduces to equation (|14D and equation (^2]) is equivalent to equation 

To generalize equation (|I0|) to polarization the Legendre tensors should be constructed 
in a different way. Polarization depends both on n and on (ei, 62) in the plane perpendicular 
to it. One may form the linear combination m = 2^^/^(ei + ^62) and construct the 
appropriate tensors by combining n and rh. We expand the polarization perturbation as 

Ap = E(2/ + 1) 2AM-P)-^C[[k)-'G^,,...,, ,Pr''. (24) 
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These "spin" Legendre tensors are symmetric combinations constructed out of n, m and 
the metric g^^ using the recurrence 



(25) 



The hierarchy begins at Z = 2 with 2P2 = 3 rrtmK To derive the Boltzmann hierarchies for 
polarization the following properties are useful and can be proven by induction 



yn«2 I 







9hi2 2^1 

Pl-l Ui, 2r-i 

Using equations (p3D, ( pS]) and p6| ) we obtain 

G\h--ilifi^ Pl^ ' ~ 







(26) 



(/ + 3)G|.,...,,, + (/ - 2)/326f^|,,..,,_, 



{21 + 1 

These relations together with 

rffi / dV (G|,,...,, 2Pr -'' 2P^'"")\^n' = 



.(27) 



(28) 



can be used to show that our choice of normalization in equations (^31) and (^Sf) coincides 
with that of flat space. 



The Boltzmann equation 

At + n'AT\i 



6 



{h + G^i,P^)+ATi 



Thomson 



Ap + Tn}Ap\i 



A 



P\Thomson 



(29) 



now becomes a hierarchy 



Atq 
Ati 

At2 

Ati 

2Api 



kAri - - 




^ [WAto - 262AT2] + ^Tl 

3 2 - 

^ [262 Ati - 363 Ats] + -7^k%a + 

[/6,AT(i-i) - (/ + l)6i+iAT(i+i) 

^ [(/ - 2)6; 2Api_i - (/ + 3)6/+i 2Api+i] - k 2AP1 ^ 



n ^ 

— At2 

10 



kAt; , / > 2 



2/ + 1 



20 



(30) 
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where 11 = At2 ~ 12 -2Api and b = — SK/k"^ = (3b2/k. The same hierarchy (equation ^0|) 
but without Thomson scattering and polarization apphes to massless neutrinos, while for 
massive neutrinos the hierarchy depends on the momentum as well (e.g. |Ma fc Bertschinger 



1995). 



Finally, the power spectra are given by 

CiT,E)i = (47r)2|/52rf/5P(/5)|A(r,^)K/3,r = ro)p 

Cci = (Anf I P'df3P{P)ATiiP,T = To)AEiiP,r = ro) 



A 



El 



(1 + 2)1 



.A 



PI, 



(31) 



with P{l3) denoting the primordial power spectrum. Equation (|3TD only applies to flat and 
open universes, whereas for the closed universe the eigenvalues of the Laplacian are discrete 
so the integral is replaced with a sum over K~^f3 = 3, 4, 5... {K~^P = 1, 2 are pure gauge 
modes, Abbott & Schaefer 1986). The usual choice for the power spectrum is 



P(/5)cx 



PiP' - K) 



(32) 



which has equal power in the curvature perturbation per logarithmic interval of k (Lyth 



& Stewart 1990, White & Bunn 1995). Note that in equation (^) one integrates over /3 
instead of k as in the fiat case. 



4. Integral Solution 



Having derived the Boltzmann equation for temperature and polarization in the 
previous section we proceed to derive the integral solution, which was the basis of the 
numerically efficient algorithm in fiat space (Paper I). An alternative, more geometrical 
derivation of the integral solution is presented in the appendix. In the fiat case the 
temperature and polarization multipoles can be written as a time integral of the product of 
a source term and a geometrical term, which is the solution of the source-less Boltzmann 
equation. The geometrical term (given in terms of spherical Bessel functions in the fiat 
case) becomes a function of two parameters in a non-fiat universe, because of the additional 
scale in the problem, the curvature of the universe. The source term can be expressed as in 
the fiat case in terms of the photon, baryon and metric perturbations. The main property 



-14- 



of this solution remains unchanged: the source term is a slowly varying function of the 
wavenumber, while the geometrical term, which oscillates much more rapidly, does not 
depend on the specific cosmological model except through its curvature. 



To obtain the integral solution it is useful to work in spherical coordinates. The 
unperturbed metric can be written as 



where the coordinate x is related to r by 

K^^/^sinK^/'^X, K >0 
r{x) = { X, K = 

'-K)~^/^sinh{-Ky/^X, K<0 



(33) 



(34) 



In these coordinates the eigenfunctions of the Laplacian (equation 0) are given by 

Gsphix) = <^'^ix)yiU^). (35) 

The radial functions ^^^(x) are the so called ultra-spherical Bessel functions. From the 
expression for the Laplacian in spherical coordinates it follows that ^^^(x) obey the following 
differential equation (Abbott & Schaefer 1986), 



d^u\ 



+ 



rixV 



u, 



0, 



(36) 



where u^(x) = ^'{x)^^i3{x)- the flat case the solutions for ^^^(x) reduce to the famihar 
spherical Bessel functions ji{kx)- The derivative of an ultra-spherical Bessel function can 
be written as (Abbott & Schaefer 1986), 



21 + 1 



ibMx) ~ {I + i)bi+iK+\x) 



(37) 



where the derivative is with respect to x = ^ t, so that dx = —dr. Comparing equations 
(0) and (^) we see that the ultra-spherical Bessel functions and their derivatives are 
solutions of the Boltzmann hierarchy in the absence of scattering and gravity. One can thus 
make the following ansatz to solve the Boltzmann hierarchy 



ATiir) = / dr'e-^^^'^'^ (r - r')5o(r') + (r - r')5i(r') + (r - r')52(r' 



(3^ 



Here coefficients Sq, Si and 5*2 were introduced because there are sources for Z = 0, 1, 2 in 
equation (RDI). By simple substitution of equation (RSI) in (150) we obtain. 
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3kU 



(39) 



The values of the ultra-spherical Bessel functions and their derivatives for x = are also 



derivatives can be obtains from equation ( P^ ) using ^^^(O) 



needed in the calculation, $^3(0) = 5/o, '^'^^(O) = |5n and $^3(0) 



%k%5n - ^Sio. These 



15 
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Finally, integrating equation (|38|) by parts one can eliminate the derivatives of 
ultra-spherical Bessel functions and the solution can be written in the following form 



ST{f3,T) 



/■TO , 

/ tifo n 3n ^ 

-K/ • "X ■ f'^b 311 ' 



+ 



(40) 



where we have defined the visibility function g = ke^'^. Comparing this expression with 
the equivalent one for the flat universe (Paper I) one can see that the only difference in the 
present case is that 11 is always divided by b, which comes from the fact that $^(0) = Y^k% 
while 72(0) = -f^k"^- Thus to transform from the flat to the general solution one should 
replace the Bessel functions with their non-fiat generalization and introduce a factor of b to 
account for the different normalization of these functions at the origin. 

To solve the polarization hierarchy we use the recursion relation (Abbott & Schaefer 
1986), 



where 



k2 _ K{1 + l)2$'+i(x) = {21 + 1) cot;,(x)$;3(x) - Vk^^KPK-\x) 



' K^l^ cot K'^/^x. K >Q 



(41) 



cot^(x) = < X-\ K = Q 



(42) 



^ {-Ky/^ coth{- Ky/^x, K<0. 



Using this and equation ( pT]) one finds 



21 + 1 



i-i 



{l-2)bi^-{l + 3)bi+i- 



i+i 



(43) 
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where the time derivative is taken on ($^/r^). This means that $^/r^ solves the polarization 
free streaming hierarchy equation ( pO] ) and is therefore the Green's function for polarization. 
The full solution is obtained by a substitution into equation (PU]): 

An alternative derivation of this result and photon transport in general can be found in the 
appendix. 

Equations (^) and (Q) are the main results of this section. The numerical 
implementation of these solutions lead to a significant reduction in computational time in 
the flat case. For this to be a successful strategy we need to show that the ultra- spherical 
Bessel functions can be computed efficiently. We turn to this subject next. 



5. The Ultra-Spherical Bessel Functions 



The main difficulty in the numerical implementation of equations pO|) and (|4^) is 
the calculation of the ultra-spherical Bessel functions The method used in the flat 

case was to precompute and tabulate these functions for all the values of interest. This 
was feasible because the functions only depend on the product kx- In the non-flat models 
another scale is introduced in the problem, the curvature of the universe. The ultra-spherical 
Bessel functions become functions of both (3 and x separately. Their tabulation is not 
feasible because these functions are rapidly oscillating in both parameters, and an excessive 
amount of memory would be needed to store the functions sufficiently densely to assure 
accurate results. In this section we develop an efficient method to compute the functions 
rapidly and with the necessary accuracy for the purpose of CMB calculations. 

The simplest approach for calculating the ultra-spherical Bessel functions would be 
to use the recursion relation in equation (^T]). This method is usually recommended as 
the most efficient if only one particular value of is required (e.g. Press et al. 1992). 

Unfortunately, using this recursive method results in computational time significantly 
longer than the time needed to compute the source term, so that the total integration time 
becomes excessive. The main drawback of this method is that for every time step needed 
to compute the time integration in equation (|40|) all are calculated up to the required 
/max- However, the smoothness of the Ci spectra allows us to calculate them sparsely and 
still obtain very accurate results. The structure of the integral solution in equation ( ^01 ) 
does not couple different / modes, and one can use this property to signiflcantly reduce 
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computational time (Paper I). The second shortcoming of the recursive method is that one 
does not use the value of $^ at the previous time step, although to guarantee accurate 
integration of equation (^Op one needs to sample $^ sufficiently densely that differences 
between successive values are small. Both shortcomings suggest that one could efficiently 
compute the ultra-spherical Bessel functions directly from the differential equation (^6l). 
Although this is not as efficient as the recursive method for any given time r, it is much 
more efficient than the recursive method if one requires the values of the ultra-spherical 
Bessel functions for a number of time steps between and tq. Moreover, one only needs to 
evaluate the functions at the needed values of /, which is typically every 50th value if 1% 
accuracy is desired. 

The structure of the ultra-spherical Bessel functions is oscillatory and similar to the 
ordinary Bessel functions (figure Q). The differential equation ( |5BD can be viewed as the 
Schrodinger equation for a particle with energy E = jS"^ in a potential V = l{l + l)/r(x)^. 
Thus for X such that /5r(x) > ^1(1 + 1), where the "energy" is greater than the "potential", 

the solution is oscillatory. For /?r(x) < ^1(1 + 1) there is a growing and a decaying solution, 

oc r^"*"^ and oc r~\ respectively. The growing solution corresponds to ^^(x)- In 
order to obtain an accurate numerical convergence the integration needs to be started in 
the regime where ^''/^{x) is small, which in our case means starting the integration close to 
/3r{x) ~ I- The equation then needs to be evolved in the direction of increasing x until 
recombination, where x ^ ^o- The integration of $^(x) therefore proceeds in the opposite 
time direction than the evolution of Boltzmann, fluid and Einstein equations (§2). If one 
were to start the integration at early times and evolve the Bessel functions to smaller radial 
distances x (i-e- the present time), the integration would be numerically unstable, as the 
decaying uj^ oc r~' mode would increasingly contaminate the solution. 

The starting value for the integration of equation ( ^6]) has to be chosen so that 
$^(x) has some flxed small value, typically 10^^. The function is changing rapidly for 
/3r(x) < + 1), so one cannot start just anywhere below this value, because one would 
soon be in the regime where $^(x) is excessively small and numerical round-off errors 
become untolerable. One could use equation (^1]) as a downward recursion for a given value 
of (3 where $j3(x) is of order 10~^, but this leads to two difficulties. One is that where 
$^(x) = 10~^ the asymptotic expansion for $^(x) is not necessarily vahd as x niight be 
too large, so starting values of x based on the asymptotic expansion do not necessarily lead 
to the desired value of the ultra-spherical Bessel function. More importantly, using the 
recursion relation only for one value of x at each (3 still results in an excessive computational 
time. The solution is to precompute the starting values for integration on a grid of f3 and 
/ and interpolate between them for any given value. Because the ultra-spherical Bessel 
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Fig. 1. — Comparison between ultraspherical Bessel function $^ (x) and spherical Bessel 
function jiQo{P sinx x) as a function of /Jsin^x. The two functions agree quahtatively, but 
oscillates much more rapidly and has a higher amplitude at first oscillation. 



functions are not oscillatory in this regime accurate interpolation can be achieved with only 
a small number of precomputed values, typically 25 values of (3 for each I. 

One can further reduce the computational time by using an asymptotic expansion for 
large radial distances x, which is obtained by the WKB approximation applied to equation 



u. 



sin[e(x)] 



where 



(45) 



(46) 



Here e is a constant phase that can be obtained analytically, but we choose it in such a way 
as to match the phase at the point where we switch from integrating equation (^) to the 



WKB approximation (equation^). Sufficiently high accuracy is achieved by switching to 
WKB approximation after one hundred oscillations of m^. 
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6. Approximations 



In this section we compare the results of exact calculations with some of the 
approximations, both those already used in the literature and some based on the integral 
solution. The goal of this section is to estimate the accuracy of flat space approximations, 
which avoid the need to compute the ultra-spherical Bessel functions and were often used in 
the context of reconstruction of cosmo logical parameters from the CMB. To avoid dealing 
with the complicated structure of ultra-spherical Bessel functions one can compute the Ci 
spectrum for a flat model and then rescale it in / as a result of the change in the angular 
diameter distance to recombination and in the size of the acoustic horizon (e.g. [Bond et 



al. 1994| , |Jungman et al. 1996| , Hu & White 1996). The quantities determining the plasma 
sound horizon and its time evolution are Q^h'^ and Vtmh'^ ■, with fif, and VL^ the density of 
baryons and matter in units of critical density. If these two parameters are kept fixed 
the only change in the spectra for small angular scales should arise from the change in 
the size of the acoustic horizon. On large angular scales this approximation fails because 
of additional effects from the decay of the gravitational potential in VLm 7^ 1 universe 
(integrated Sachs- Wolfe contribution), as well as curvature effects on the initial spectrum 
and on the radial eigenmodes. The simplest and most often used approximation is to 
rescale the spectrum by VLq^^"^ (e.g. Jungman et al. 1996), while keeping Vt^ah'^ and ^b/i^ 
fixed. The dependence of the acoustic horizon on VLq is only poorly approximated by this 
scaling and leads to significant errors, as shown in figure ^ for the case VLq = 0.2. Clearly 
this approximation is too inaccurate to be useful for numerical analysis. 

Significantly better is the approximation where the angular diameter distance to the 
last-scattering surface r{xr) is correctly calculated and then used to rescale the spectrum 
with (Hu & White 1996), 

Xr 

where Xr is the radial distance to the recombination epoch. Recombination occurs roughly 
around z ~ 1100, but there is some uncertainty in its exact value and it varies somewhat 
with fifo. The normalization of the spectrum is also a problem, since this scaling is only 
valid for high enough values of / making it impossible to use the lower multipoles for 
normalization. To avoid the two problems mentioned above we performed the full open 
calculation and then shifted the fiat spectra so that they agreed exactly at the first peak. 
This is therefore the most favorable case, and the above mentioned difficulties will make 
the agreement worse. Figure |^ shows this approximation together with the exact results. 
It obviously fails on large scales because of the effects mentioned above. While some of 
these effects can be modeled analytically it is not obvious how one connects the small 
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Fig. 2. — Comparison between the full calculation (solid line), fig rescaling of the flat 
spectrum (dotted line) and the angular diameter rescaling of the flat spectrum (dashed hne) 
for the case with Qq = 0.2. The latter approximation fares significantly better than the 
first one on small scales, while on large scales additional effects in open universes make the 
agreement poor for both approximations. 



scale and large scale approximations. On small scales the agreement is significantly better 
and the approximation deviates from the exact calculation by about 3-4%. While this 
approximation fares significantly better than the other approximation above, it is still not 
sufficiently accurate for the forthcoming high precision data. 

We also tried a more sophisticated version of the above approximation which solves 
some of the problems. In this approximation one uses the integral solution and calculates the 
source term using equations presented in §2, but the radial eigenfunctions are replaced with 
their fiat space approximations. This means replacing the ultra-spherical Bessel functions 
in equation (BO) by the spherical Bessel functions with the argument kr{x), so that the 
appropriate angular diameter distance relation is used. The small scale normalization and 
epoch of recombination are correctly calculated in this case and the integrated Sachs- Wolfe 
term is included. If this were a good approximation one could reduce the computational 
time by being able to tabulate the geometrical term as was done for the fiat models (Paper 
I). While replacing the ultra-spherical Bessel functions with the fiat approximation gives 
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qualitative agreement with the exact results, it does not converge to the exact solution in 
the limit of large /. This is shown in figure |I|, where the two are compared for / = 100. 
The approximation gives correct qualitative behavior and properly describes the position 
of the transition from ~ to the oscillatory regime, as well as the amplitude of 

oscillations. If the last-scattering surface is thin, then in the time integration of equation 
(^) the spherical Bessel function does not change and it can be taken out of the integral. 
Then in equation (|3T|) the square of the Bessel function is multiplied by the square of the 
source averaged over time, this product is integrated over the primordial spectrum. Because 
both the primordial spectrum and the source are slowly changing with k, this essentially 
amounts to an averaging of jf over its oscillation period. The result only depends on the 
amplitude and not on the oscillation frequency of the spherical Bessel function. This is the 
underlying reason why the approximation of rescaling the distance to the last-scattering 
surface gives qualitative agreement in the final spectrum and in particular, it correctly 
predicts the positions of the acoustic oscillation peaks. 

However, as shown in figure |l| the frequency of oscillation is different for the two 
functions, and this has important consequences when the finite width of the last scattering 
surface is taken into account. In this case one has to explicitly integrate the sources over 
the recombination epoch as in equation (^Ol). The integrated Sachs- Wolfe term cannot be 
calculated correctly using this approximation because the time dependence of the Bessel 
functions is important. This disagreement is most severe for the lower multipoles. The 
agreement for the part proportional to the visibility function is not good either, because 
ji[(3 sinh(x)] oscillates faster in x than (the phase is proportional to (3 sinh(x) rather than 
to Px), so the damping is more severe than in its ultra-spherical counterpart. The situation 
is even worse for the terms proportional to the first and second derivatives of the visibility 
function (equation ^0[). In the limit where the visibility function can be approximated by 
a 5-function this is easy to understand. The derivatives of the visibility function act as 
derivatives of the (5-function so these terms involve derivatives of the ultra-spherical Bessel 
functions, which are not well approximated by the derivatives of ji{(3 sinhx)- This can be 
seen from the fact that the latter, although having the same amplitude as $^3, oscillate 
more rapidly. This problem could be solved by using equation (37) to express in terms 
of the and ^''^^ and approximating these by their fiat counterparts. But the rest of 
the problems mentioned above will remain. This fiat space approximation is therefore not 
suitable for the integral approach and its results are even worse than the simple rescaling of 
the spectrum with the angular diameter distance. To summarize the discussion, none of the 
fiat space approximations can match the exact results at better than 5% accuracy, which is 
not sufficient for the expected future observational precision. 
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7. Implementation 

We are now in a position to present an efficient implementation of the method 
developed in the previous sections. Some of the technical issues of this implementation 
are identical to those in the fiat case (Paper I) and will not be repeated here, but non-flat 
geometry also introduces some additional complications. The temperature spectra are 
calculated by integrating equation P0|), while the source in this equation is obtained by 
integrating the equations presented in §2. One of the advantages of calculating the Ci 
spectra this way is that the source in the integral only depends on the ffist few photon 
multipoles, so accurate results can be obtained truncating the Boltzmann hierarchy at a 
very low value of /. As discussed in Paper I this is only possible if one uses a non-refiective 
boundary conditions to close the hierarchy in equation (pOD, so that the propagation of 
power from large scales to small scales is not reflected back onto the lowest moments of 
photon distribution. In the absence of scattering and integrated Sachs-Wolfe contribution 
equation (^) becomes identical to equation (|37D, hence one can use the recursion relation 
in equation ( ^T] ) to relate A^+i to A; and thus close the hierarchy. In reality this relation is 
not exact, but since only the lowest moments need to be accurately calculated, one can use 
this closure scheme at some moment /q, which is sufficiently high so that any reflection of 
power from Zq down to the lowest moments will be negligible. In practice this is achieved 
to a high accuracy for /q = 7, similar to the flat case. This therefore reduces the number of 
equations needed to be solved from a few thousand in the standard Boltzmann code to a 
few dozen in the integral approach, making the calculations significantly faster. 

It was mentioned in the previous section that the sources and the ultra-spherical Bessel 
functions are calculated in the opposite directions of time, hence the integral in equation 
(1^) cannot be computed at the same time as the sources are computed. The radial 
eigenfunctions in equation (|^) are calculated from equation ( ^6]) using a finite differencing 
scheme such as Runge-Kutta 4th order integration, a procedure that automatically produces 
sufficient sampling points in time for an accurate integration. Because the sources do not 
depend on / and as discussed below need to be computed at far fewer points than the 
radial eigenfunctions, we compute them first at selected values of wavevector k by solving 
the system of equations presented in §2. Their values are stored at selected time intervals. 
During the integration of equation ( PBD the sources are interpolated using cubic splines from 
their tabulated values and equation (jl^) is solved simultaneously. To obtain an accurate 
interpolation of the sources, at most a few hundred values in time are needed, covering the 
time interval from recombination until today. 

The main advantage of the integral method is in the sampling of the wavenumber k for 
the source term. The sources vary typically on a scale ~ l/r,. where is the conformal 
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time of recombination. The oscillations in are approximately 50 times faster, as they 
occur on a scale of the order of the inverse of the present conformal time. ro. We take 
advantage of this by calculating the sources at fewer values of k than what is needed to 
properly sample the variations of the radial eigenfunctions, and then interpolate them to 
obtain their values for the rest of the values of k. This reduces the number of k modes 
where the system of equations in §2 is solved, and is an important factor that enhances the 
performance of the method. 

Because the Ci spectra are so smooth, they can be sampled in every 50th /. except 
for low values of / where denser sampling is needed. The rest of the C/'s arc obtained by 
interpolation. This reduces the number of radial eigenfunctions that need to be computed, 
so that approximately 45 values of I are needed up to l^ax — 1500. They need to be sampled 
at least at a few points per oscillation, resulting in approximately 2000 points for the above 
Imax- One therefore needs to solve 10^ differential equations for the radial eigenmodes 
and only about 3000 differential equations for the source term. However, the differential 
equations for the sources are more complicated, so the actual computational time for the 
two terms is comparable. Typical CPU time for an open model with 1% accuracy up to 
^max = 1500 is about 60 seconds on an ALPHA-500 workstation, which is about two times 
slower than for a comparable model in a fiat universe and orders of magnitude faster than 
any other method. 



8. Discussion 

We derived the photon Boltzmann hierarchy for scalar temperature and polarization 
perturbations in a general Robertson- Walker universe using a new expansion method in 
spin-s harmonics. Recently, a similar expansion has been presented by Hu & White (1997) 
and applied to scalar, vector and tensor harmonics, but only in a fiat universe. Combining 
these techniques will provide a complete treatment of cosmological perturbation theory in 
a perturbed Robertson- Walker universe and will be presented elsewhere (Hu et al. 1997). 
In this paper we concentrate on providing an efficient method to compute scalar CMB 
temperature and polarization anisotropics in a non-flat universe, generalizing the methods 
presented in Paper I. For this purpose we present an integral solution to the Boltzmann 
hierarchy equations. Temperature and polarization perturbations are written using Green's 
method as a time integral over the source term multiplied by a radial eigenf unction. The 
first term depends on the cosmological model, but not on the multipole moment, while 
the latter depends only on the curvature of space, but not on the rest of the cosmological 
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parameters. This split clearly separates between effects that depend on the source and those 
that depend on geometry and is more physically intuitive than the differential formulation. 
The numerical implementation of the integral method requires an efficient and accurate 
way of calculating the ultra spherical Bessel functions, because these cannot be stored in 
advance as in the fiat case. We achieve this by solving their differential equation. The 
implementation of the integral solution leads to a fast and accurate method for computing 
CMB anisotropy and polarization spectra with computational times of the same order as in 
the flat case. It should be noted that once the ultra spherical Bessel functions are calculated 
computing vector and tensor spectra poses no further calculational difficulties. 

The integral method presented in this paper can be used for studying theoretical 
predictions from various cosmological models with the purpose of comparing them to the 
real data. The implementation of the method is pubhcly available. The output of the 
calculation consists of both the temperature and polarization spectra as well as the matter 
density power spectra, so one can use it to analyze both clustering and CMB data. In 
addition, an accurate small scale normalization (erg) to the COBE data is provided by using 
the 4 year analysis of the power spectrum (Bunn & White 1997). As such it should replace 
various approximative methods used in the literature if high precision is required. A new 
generation of experiments will be able to determine both the density power spectrum (e.g. 
Sloan Digital Sky Survey and the 2dF survey) and the CMB spectrum (MAP, Planck) to 
an exquisite accuracy. The integral method presented here provides a tool for rapid and 
highly accurate analysis of theoretical models over a large range of parameter space. This 
should be useful in achieving a higher level of precision in the determination of the true 
cosmological model. 

We acknowledge useful discussions with W. Hu, M. Mahachek and M. White. We 
also want to thank the referee, Douglas Scott for many useful comments. This work was 
partially supported by grant NASA NHG5-2816. 



A. Appendix 

In this appendix we present an alternative and more intuitive derivation of the integral 
solution. We first present a solution of the photon free streaming equation and then 
construct the general solution from it using Green's method. This method clarifies the 
structure of the radial part of the integral solution and the connection between spatially 
varying perturbations and the resulting angular variation in the sky. 
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Polarization is described hjQ±iU which have spin s = ±2. The fields T, Q ±iU 
are function of position in space specified in spherical coordinates by (x, 0, 0) as well 
as direction in the sky. We take the observer to be located at the origin, in spherical 
coordinates at x = so we are interested in T and Q ±iU a.s a function of the direction of 
observation n at the origin today. As discussed in the main text polarization is not only a 
function of n but also of the two directions perpendicular to n, (61,62), which are used to 
define the Stokes parameters. Thus, 

T = T{x,0,(f),h,T) 
Q±iU = (g±zt/)(x,^,0,n,ei,e2,r) (Al) 

We can make a single vector out of these quantities X = {T,Q + iU, Q — iU). 

We begin by considering the solution to the free streaming equation, 

and then construct the complete solution from that. Notice that this is only true if 
the polarization vector does not rotate relative to {eg, e^i)) during propagation, which 
can happen in some Bianchi models ( [Iblman fc Matzner 198^ ). Given the field 



9, 0, n, ci, 62, T*) = X* at some initial time r* the free streaming solution for the field 
at the origin today is simply 

X{Q, h, ee, e^, tq) = X*{x*, 6', 0, n, eg, e^, r*) (A3) 

where {9, (p) specify the direction of n in spherical coordinates and we use {eg, e^) = (ei, 62). 
Because photon position and time are related via x* = tq — t* we may drop the explicit 
dependence on r* in the following. When considering the perturbations produced by one 
single scalar mode, G{x), the complete initial field X* = {T*, Q* + iU* , Q* — iU*) must be 
constructed from that scalar mode. The temperature is a spin zero function of n and so 
G and its radial derivatives can be used to construct the most general form of X* for this 
mode, 

AT* = oeoG + oeiGi.ri* + . . . (A4) 

where we introduced the coefficients to expand the temperature perturbation in terms 
of the mode functions and their derivatives. We concentrate on the first term, which 
corresponds to an initial distribution that is spatially varying as G and locally isotropic. It 
is through the free streaming of the photons that the spatial variations of the source create 
an angular dependent distribution function despite being isotropic initially. Taking G as 
the eigenfunction of the Laplacian in spherical coordinates Ggpher = ^''i3{x)'^im{9 , 0) ( [Abbott 
fe Schaefier 1986| ) we find using equation ([A3|) , 



AT(0, n. To) = oeo$^(ro - r*)YUO, 0) (A5) 
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This means that $^(to — r*) is a solution of the free streaming equation, which we used 
in the main text to construct the full solution of the Boltzmann equation for temperature. 
The source terms for temperature consist not only of isotropic term, but also of higher 
multipoles. Therefore, for the complete solution we have to include also the terms 
proportional to the derivatives of the ultra-spherical Bessel functions, corresponding to 
higher order terms in the expansion of equation 



For polarization we need an initial field of spin 2, which is a function of n and 
{eg,e^). We can construct a function like this by contracting the derivatives of G with 
m = 2~^/^(ee + ie^), fn = 2~^/^(ee — ie^) and n in the following way, 

A{Q + ill)* = 2^0 Giijrh'rh^ + 2^iG\ijkm'm^h^ + . . . 

A(Q - lUy = _2eo G\ijm'rfi' + ^2tiG\ijkmrni h'' + ... (A6) 

where again we have introduced expansion coefficients ±2^1- For polarization we only need 
to consider the first term in each expansion because this suffices for the angular structure 
of the source. The free streaming solution is 

A((5 + if/) (0,n, 60,60) = G\ijrh}rn> 

A{Q - iU){O,n,e0,e^) = -2eo G\ijmm^ . (A7) 

By taking derivatives of Gspheri^) we can show that 

• • I $Uy) 

Gspher\^JM^ = ^ {I + 2)\ / {I - 2)1-^ 2YUd A) 

G,pHer\^Jmrh' = + 2) ! / (/ - 2) ! -g^ _2>^„ (e, 0) (A8) 
so that equation ([A7|) becomes 



/ $Uy) 

A(Q + ^f/)(O,n,ee,e0) = 2eo J (l + 2)1/(1 - 2)\^^ 2Y1UO, <!>) 

A{Q-tU){0,h,ee,e^) = _2eov/(/ + 2)!/(/ - 2)!-g^ _2r,™(^, c^), (A9) 

with X = To ^ T- This shows that $^/r(x)^ is the free streaming solution for polarization as 
was obtained using a different method in the main text. 

The complete Boltzmann equation for polarization is of the form 

i^i^El = _,a(q->c/)-|5g,,™a' (Aio) 
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The G\ijm^m^ leads to the usual (1 — /i^) oc (1 — -P2(/u)) dependence of Thomson scattering. 
The angular dependence of the source is the simplest spin two function that can be 
constructed from a scalar. One can interpret the above equation as stating that at each 
time r* a contribution d{^Q + zAf/) oc kIlG\ijm^m^dT is generated and then free streamed 
until today. If there is scattering along the way only a fraction exp(— k) of the photons 
reach the observer. Given the free streaming solution (equation \K3i ) the complete solution 
can be written as 



A{Q + tU)iO,n,ee,e^,To) = (I + 2)\ / {I - 2)\ 
A(g-«f/)(0,n,ee,e<^,ro) = -J {I + 2)\/{l - 2)\ 



drg{T 



3n<l>;3(ro 



(Aii: 



where gir) is the visibility function. When considering solutions of equation ( pOD 
corresponding to one mode that is the generalization of a single plane wave a superposition 
of solutions ( All ) for all / with m = must be used. 



The above solutions for temperature and polarization help to illustrate the relation 
between the radial part of the mode function G and the free streaming solution for A" as a 
function of angle in the sky. When we look in a given direction on the sky we observe what 
was happening away from us at a distance X = tq — r*. This couples the spatial variation of 
the mode functions with the angular variations of the observed distribution. The angular 
dependence of the sources, which is constrained by the spin nature of the variables, is also 
responsible for the specific form of the integral solution. The source for scalar polarization 
has to be constructed out of derivatives of the mode functions which change the radial 
dependence of the solution from to $^/r^. 

These arguments can be extended to vector and tensor modes. The radial and angular 
dependence of the mode functions differs from the scalar case and can be found in Tomita 
(1982). In the case of temperature the source is proportional to /i^r = hijfi^n^ , while for 
polarization it is proportional to hijrftm^ and h 



- * - J 
m m 



for Q + ill and Q — ill for any type 
of perturbations. This together with the explicit mode functions (Tomita 1982) can be used 
to verify the fiat space solution obtained previously (Zaldarriaga & Seljak 1997). Detailed 
analysis of vector and tensor non-flat Boltzmann equation will be presented elsewhere (Hu 
et al. 1997). 



REFERENCES 



Abbott, L. F., & Schaefer, R. K. 1986, ApJ, 308, 546 



- 28 - 



Bardeen, J. M. 1980, Phys.Rev.D, 22, 1882 

Bertschinger, E. 1996, in "Cosmology and Large Scale Structure", ed. R Scliaffer et. al. 
(Elsevier Science, Netherlands), p. 273 

Bond, J. R., et al. 1994, Phys.Rev.Lett, 72, 13 

Bond, J. R., Efstathiou, G., & Tegmark, M. 1997, |astro-ph/9702T00| 
Bond, J. R., & Efstathiou, G. 1984, ApJ285, L45 

Bucher, M., Goldhaber, A. S., & Turok, N. 1995, Phys.Rev.D, 52, 3314 

Bunn, E. P., & White, M. 1997, ApJ, 480, 6 

Chandrasekar, S. 1960, Radiative Transfer (Dover, New York) 

Jungman, G., Kamionkowski, M., Kosowsky, A., & Spergel, D. N. 1996, Phys.Rev.D, 54, 
1332 

Goldberg, J. N., et al. 1967, J. Math. Phys. 8, 2155 

Gorski, K. M., et al. 1965, ApJ446, L67 

Gouda, N. & Sugiyama, N. 1992, ApJ, 395, L59 

Hu, W., & White, M. 1996, ApJ, 471, 30 

Hu, W., Bunn, E. & Sugiyama, N. 1995, ApJ447, L59 

Hu, W. and White, M. 1997, |astro-ph 9702T70 



Hu, W., Seljak, U., White, M., and Zaldarriaga, M. 1997, in preparation 
Kamionkowski, M., Spergel, D. & Sugiyama. N. 1994, ApJ, 426, L57 
Kamionkowski, M., Kosowsky A., and Stebbins A. 1997, Phys.Rev.D, 55, 7368 
Kochanek, C. S. 1996, ApJ, 473, 595 
Liddle, A. R. et al. 1995, M.N.R.A.S., 278, 644 
Lifshitz, E. M. 1946, J. Phys. USSR, 10, 16 
Lineweaver, C. H., Barbosa, D., Blanchard, A., and Bartlett, J. G. 1996, |astro-ph /9610133 
Lyth, D. H. & Stewart, E. D. 1990, Phys. Lett., B252, 336 



- 29 - 



Ma, C.-R, & Bertschinger, E. 1995, ApJ, 455, 7 

Peebles, P. J. E., 1993, "Principles of Physical Cosmology", Princeton University Press, 
Princeton 

Peebles, P. J. E., & Yu, J. T. 1970, ApJ, 162, 815 



Perlmutter, S., et al. 1996, preprint |astro-ph/9608192" 



Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical 
Recipes in Fortran, 2nd ed. (New York: Cambridge University Press) 

Ratra, B. & Peebles P. J. E. 1994, ApJ, 432, L5 

Rocha, C, and Hancock, S. 1996, ^stro-ph/96lT228 



Seljak, U., & Zaldarriaga, M. 1996 (Paper I), ApJ, 469, 437 
Smoot, G. et al. 1992, ApJ, 396, LIS 
Tomita, K., 1982 Prog. Tlieo. Phys., 68, 310 

Tolman, B. W. and Matzner, R. A. 1984, Proc. R. Soc. Lond. A 392, 391 

White, M., & Scott, D. 1996, ApJ, 459, 415 

White, M., & Bunn, E. 1996, ApJ, 460, 1071 

Wilson, M. L., & Silk, J. 1981, ApJ, 243, 14 

Zaldarriaga, M., & Seljak, U., 1997 Phys. Rev. D 55, 1830 

Zaldarriaga, M., Spergel D. & Seljak, U., |astro-ph/ 9702157 



This preprint was prepared with the AAS M5t;X macros v3.0. 



